—- message=FALSE——————————————————

#setwd(“D:/ET_Hands_On”)

library(raster)
## Loading required package: sp
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-27, (SVN revision 1148)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/r-library/3_6_3/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/r-library/3_6_3/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/r-library/3_6_3/rgdal/proj
library(ETUAS)
R <- raster("RETAU2_PO_08282018_50_B6.tif")
plot(R)

aoi <- createAoi(topleft = c(344200, 5184000), bottomright = c(344270, 5183880), EPSG = 32611)
r <- crop(R,aoi)

————————————————————————

csvfile <- "Othello_8_28_18_weather.csv"
MTLfile <- "08_28_18_50m_MTL.txt"
#### elev in 'm' MSL ###
WeatherStation <- read.WSdata(WSdata = csvfile, date.format = "%d/%m/%Y", lat= 46.03, long= -119.53, elev=269.443, height= 2.2,
                              columns=c("date" = 1, "time" = 2, "radiation" = 3,"wind" = 4, "RH" = 6, "temp" = 7, "rain" = 8), MTL = MTLfile)
## Warning in read.WSdata(WSdata = csvfile, date.format = "%d/%m/%Y", lat =
## 46.03, : As tz = "", assuming the weather station time zone is America/
## Los_Angeles
image.DN <-loadRawSpectra(aoi=aoi)
## Warning in `[<-`(`*tmp*`, i, value = raster(list.files(path = path, pattern =
## paste0(image_pattern, : implicit list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = raster(list.files(path = path, pattern =
## paste0(image_pattern, : implicit list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = raster(list.files(path = path, pattern =
## paste0(image_pattern, : implicit list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = raster(list.files(path = path, pattern =
## paste0(image_pattern, : implicit list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = raster(list.files(path = path, pattern =
## paste0(image_pattern, : implicit list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = raster(list.files(path = path, pattern =
## paste0(image_pattern, : implicit list embedding of S4 objects is deprecated
plot(image.DN)

DEM<-prepareDEM(extent=image.DN)
## Warning in projectRaster(SRTMmosaic, destino): input and ouput crs are the same
## ---- surface properties ##------------------------------------------------------
surface.model <-METRICtopo(DEM)
surface.model$Slope[is.na(surface.model$Slope[])] <- 0
surface.model$Aspect[is.na(surface.model$Aspect[])] <- 0

## solar angles ##
solar.angles.r <- solarAngles(surface.model = surface.model,WeatherStation = WeatherStation, MTL = MTLfile)

## Incoming shortwave radiation ##
Rs.inc <- incSWradiation(surface.model = surface.model, solar.angles = solar.angles.r, WeatherStation = WeatherStation)
plot(Rs.inc)

#writeRaster(Rs.inc, filename="ISWR_PO_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)

## Surface albedo ##
albedo <- (image.DN$B)*0.237 + (image.DN$G)*0.239 + (image.DN$R)*0.118 + (image.DN$NIR)*0.098 + (image.DN$RE)*0.306
#plot(albedo)
#writeRaster(albedo, filename="albedo_PO_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)
## ---- LAI default method---- #
LAI <- LAI(method = "metric2010", image = image.DN, L=0.1)
plot(LAI)

#writeRaster(LAI, filename="LAI_PO_METRIC_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)
## ---- surface temperature with calibration and blackbody emissivity correction------#
Ts <- (1.164*(image.DN$Thermal)-7.204)*((1/0.95)^0.25) +273.15
plot(Ts)

#writeRaster(Ts, filename="Ts_PO_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)

## Outgoing longwave radiation from surface ##
Rl.out <- outLWradiation(LAI = LAI, Ts=Ts)
plot(Rl.out)

#writeRaster(Rl.out, filename="OLWR_PO_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)

## Incoming longwave radiation from atmosphere ##
Rl.inc <- incLWradiation(WeatherStation,DEM = surface.model$DEM,solar.angles = solar.angles.r, Ts= Ts)
plot(Rl.inc)

#writeRaster(Rl.inc, filename="ILWR_PO_METRIC_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)

##--- Net radiation--- ##
Rn <- netRadiation(LAI, albedo, Rs.inc, Rl.inc, Rl.out)
plot(Rn)

#writeRaster(Rn, filename="Rn_PO_UAS_METRIC_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)
## ----Soil Heat Flux ---- ##
G <- soilHeatFlux(image = image.DN, Ts=Ts,albedo=albedo,Rn=Rn, LAI=LAI)
plot(G)

#writeRaster(G, filename="G_PO_UAS_METRIC_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)
## ----Momentum roughness length ----- ##
Z.om <- momentumRoughnessLength(LAI=LAI, mountainous = TRUE, method = "short.crops",surface.model = surface.model)
plot(Z.om)

#writeRaster(Z.om, filename="Zom_PO_UAS_METRIC_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)
## Hot and cold Anchor pixels ##
hot.and.cold <- calcAnchors(image = image.DN, Ts = Ts, LAI = LAI, plots = F,albedo = albedo, Z.om = Z.om, n = 5, anchors.method = "flexible", deltaTemp = 5, 
                            verbose = FALSE)
## Warning in calcAnchors(image = image.DN, Ts = Ts, LAI = LAI, plots = F, : I can
## only find 1 anchors with cold pixel conditions
## Warning in calcAnchors(image = image.DN, Ts = Ts, LAI = LAI, plots = F, : I can
## only find 1 anchors with hot pixel conditions
#write.csv(hot.and.cold, file = "Anchors_PO_UAS_METRIC_8_28_18_50m.csv")


## Sensible heat flux ##
H <- calcH(anchors = hot.and.cold, Ts = Ts, Z.om = Z.om, WeatherStation = WeatherStation, ETp.coef = 1.05,
           Z.om.ws = 0.03, DEM = DEM, Rn = Rn, G = G, verbose = FALSE)
## Warning in calcH(anchors = hot.and.cold, Ts = Ts, Z.om = Z.om, WeatherStation =
## WeatherStation, : u200 less than threshold value = 0m/s. using u200 = 4m/s

#plot(H$H)
#writeRaster(H$H, filename="H_PO_UAS_METRIC_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)
#plot(H$dT)
#writeRaster(H$dT, filename="dT_PO_UAS_METRIC_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)
## Reference ET for alfalfa crop##
ET_WS <- dailyET(WeatherStation = WeatherStation, MTL = MTLfile)

## Daily ET, METRIC based ##
ET.24 <- ET24h(Rn, G, H$H, Ts, WeatherStation = WeatherStation, ETr.daily=ET_WS)

#writeRaster(ET.24, filename="Daily_ET_PO_UAS_METRIC_8_28_18_50m.tif", format="GTiff", overwrite=TRUE)